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I.  Introduction 


The  optimal  trajectory  of  a  maneuvering  hypersonic  reentry  vehicle,  such  as  the  Space 
Shuttle,  is  generated  as  the  solution  to  a  constrained  two  point  boundary  value  problem. 
The  boundary  value  problem  is  formulated  via  identifying  both  initial  and  final  conditions, 
such  as  launch  and  landing  geographic  locations.  Furthermore,  the  problem  is  constrained 
by  limiting  angle  of  attack,  aerodynamic  heating,  and  setting  upper  and  lower  boundaries 
on  altitude.  These  conditions  define  the  Hessian  from  which  the  optimal  solution  of  the 
boundary  value  problem  is  solved. 

While  analytical  methods  exist  to  solve  a  range  of  elementary  boundary  value  problems, 
the  complexity  of  problems  such  as  a  hypersonic  vehicle  reentry  trajectory  are  better  suited 
to  numerical  methods.  This  approach  is  well  established  in  dedicated  trajectory  generation 
codes  such  as  NASA-Langley’s  Program  to  Optimize  Simulated  Trajectories  (POST)  and 
NASA-Glenn’s  Optimal  Trajectories  by  Implicit  Simulations  (OTIS).  While  the  approaches 
to  solving  the  boundary  value  problem  differ,  both  of  the  above  mentioned  codes  have  been 
shown  by  Nelson4  to  produce  nearly  identical  numerical  results  in  the  case  of  minimally 
constrained  trajectories.  In  the  case  of  a  more  heavily  constrained  problem,  however,  the 
solutions  from  the  two  codes  may  begin  to  diverge  in  solution,  or  have  the  case  when  one 
fails  to  converge  on  a  solution  altogether. 

The  necessity  to  increase  the  number  of  constraints  on  a  reentry  trajectory  is  reliant  on 
the  definition  of  the  mission  of  the  reentry  vehicle.  In  the  case  of  the  Space  Shuttle,  the 
solution  space  is  constrained  by  the  physical  constraints  mentioned  above,  but  with  few  if 
any  intermediate  waypoint  constraints.  In  the  case  of  a  more  demanding  mission  profile 
additional  constraints  are  levied  on  the  problem.  These  can  include  physical  constraints 
such  as  maximizing  the  lift  to  drag  ratio  in  order  to  maximize  the  down  range  capability 
of  the  vehicle.  Additional  constraints  can  also  be  geographic  in  nature,  such  as  avoiding 
overflight  of  specified  airspace.  Due  to  these  additional  constraints  the  ability  to  converge 
to  an  optimal  solution  becomes  dependent  on  the  numeric  method  used. 

The  objective  of  this  paper  is  to  demonstrate  the  ability  of  pseudospectral  method  to 
converge  to  a  constrained  optimal  reentry  trajectory  in  ever  increasing  levels  of  constraint 
complexity.  The  generated  numeric  solutions  and  time  to  convergence  are  evaluated  to 
those  generated  by  POST  for  the  purpose  of  comparison  to  an  established  solution  method. 
DIDO5-9  also  implements  the  pseudospectral  method;  however,  GPOCS’s  implementation  of 
phases  is  more  compatible  with  interior  point  constraints,  such  as  waypoints. 
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II.  Optimization  Methods 


As  stated  in  the  introduction,  several  numerical  techniques  exist  to  solve  constrained  two 
point  boundary  value  problems  cast  as  hypersonic  vehicle  trajectory  problems.  This  section 
will  briefly  introduce  the  numeric  optimization  techniques  used  in  both  POST  and  GPOCS. 

II. A.  Optimization  in  POST 

The  optimization  method  available  in  POST  is  the  accelerated  projected  gradient  method 
(PGA).  The  PGA  algorithm  is  an  iterative  technique  designed  to  solve  a  general  class  of 
nonlinear  programming  problems.  PGA  uses  a  cost  function  and  constraint  gradient  in¬ 
formation  in  place  of  a  multidimensional  optimization  problem  by  a  set  of  equivalent  one 
dimensional  searches.  With  this  approach,  the  initial  function  of  the  PGA  algorithm  is  to 
ensure  constraint  satisfaction,  however  the  terminal  phase  of  the  optimization  primarily  en¬ 
sures  a  reduction  in  the  cost  function.  As  with  the  case  of  many  optimization  approaches 
that  rely  on  gradient  information  to  arrive  at  a  global  minimum,  shallow  gradients  can  lead 
to  long  convergence  times  or  fail  to  converge  altogether.  One  additional  problem  with  this 
optimization  method  is  that  all  constraints  must  be  considered  simultaneously  during  the 
solution,  which  can  lead  to  non-convergence  in  the  case  of  a  highly  constrained  problem. 

II. B.  Optimization  in  GPOCS 

GPOCS  uses  the  pseudospectral  optimization  method  to  converge  to  a  minimal  cost  of  a 
two-point  boundary  value  problem.  As  stated  by  Jorris, 10,11  the  purpose  of  pseudospectral 
methods  is  to  approximate  the  continuous  solution  to  a  set  of  differential  equations  using 
polynomial  interpolation  through  discrete  points  or  nodes.  The  motivation  is  to  avoid  se¬ 
quential  integration  which  can  lead  to  divergence  and  may  prohibit  determining  a  solution. 
Lagrange,  Legendre  or  Chebyshev  polynomials,  shown  in  Figure  1,  may  be  used  to  satisfy  the 
differential  equations  at  a  discrete  number  of  nodes  N.  Computing  the  solution  at  the  nodes 
is  also  termed  collocation.  These  collocation  techniques  satisfy  all  the  nodes  simultaneously, 
thus  avoiding  the  pitfalls  of  integration,  especially  the  forward  and  backward  integration 
within  the  shooting  method.  To  ensure  computational  accuracy  of  the  solution  derived  us¬ 
ing  this  method,  the  spacing  of  the  nodes  is  of  critical  importance.  While  equally  spaced 
nodes  are  the  simplest  approach,  the  polynomial  interpolation  will  lead  to  large  errors  as  the 
number  of  nodes  is  increased.  To  avoid  this  phenomenon,  known  as  Runge  Phenomenon , 
Chebyshev  point  placement  as  defined  in  equation  1  may  be  used  to  distribute  the  nodes. 
The  roots  of  the  Legendre  polynomials,  or  Gauss  points,  may  also  be  used.  The  use  of  these 
non-uniformly  spaced  points  allows  for  increased  polynomial  interpolation  accuracy  as  the 
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Figure  1.  Comparison  of  Legendre  and  Chebyshev  Polynomials. 


number  of  nodes  increases. 

Xj  =  cos(^),  j  =  0, (1) 

For  example,  the  benefit  of  using  Chebyshev  node  spacing  over  equal  spacing  is  depicted 
in  Figure  2.  One  area  of  concern  is  that  the  midpoint  regions  along  a  trajectory  can  have 
relatively  low  node  density,  however  this  can  be  addressed  through  the  use  of  phases.  The 
term  phases  refers  to  inserting  intermediate  events  into  the  entire  span  of  nodes.  These 
additional  nodes  create  interior  boundary  conditions  by  creating  new  internal  endpoints. 
This  approach  has  the  benefit  of  adding  a  higher  density  of  nodes  in  regions  that  may 
otherwise  have  relatively  low  node  density,  which  increases  solution  accuracy.  In  practice, 
this  ability  to  segment  the  solution  of  a  set  of  differential  equations  into  phases  allows 
greater  freedom  in  defining  a  highly  constrained  boundary  value  problem.  A  more  detailed 
explanation  of  this  optimization  technique  is  covered  in  Appendix  B. 

III.  Methodology 

This  section  briefly  describes  the  methodology  used  in  this  paper  to  derive  a  technique 
that  converges  to  a  merit  optimal  solution.  The  derived  method  has  been  set  up  to  ensure 
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Figure  2.  Demonstration  of  Runge  Phenomenon  with  Equally  Space  Nodes.10 
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satisfaction  of  waypoint  and  no-fly  zone  constraints  while  requiring  minimal  computational 
time.  The  trajectory  is  constrained  by  the  equations  of  motion,  intermediate  waypoint  and 
no-fly  zone  constraints,  and  bank  angle  and  lift  control  limitations.  As  a  note,  the  bank 
angle  control  may  be  limited  by  structural  loads,  controllability,  or  heating  constraints. 
The  limitation  on  bank  angle  translates  to  a  minimum  turn  radius  to  transition  from  one 
waypoint  or  no-fly  zone  to  the  next,  while  the  coefficient  of  lift  is  used  to  modulate  the  rate 
of  descent.  Two  dimensional  analysis10  demonstrated  the  use  of  a  geometric  solution  as  an 
initial  guess  to  the  dynamic  optimization  technique.  The  work  herein  expands  upon  this 
previous  research11  to  include  a  full  three  dimensional  environment  incorporating  the  earth’s 
rotation.  This  expansion  from  a  2  dimensional  to  3  dimensional  environment  allows  for  a 
more  realistic  representation  of  hypersonic  vehicle  equations  of  motion. 

III. A.  Assumptions 

The  following  assumptions  are  made  to  scope  the  research  effort: 

1.  The  waypoints  are  specified  in  the  desired  sequence. 

2.  Inner-loop  control  is  available.  Only  the  outer-loop  or  trajectory  generation  is  ad¬ 
dressed. 

3.  Waypoints  are  sufficiently  spaced  such  that  no  two  are  within  the  turn  radius  of  the 
vehicle. 

4.  The  no-fly  zones  are  specified  as  elliptical  exclusion  zones,  radiating  from  the  center  of 
the  Earth. 

III.B.  Analysis 

The  dimensional  equations  of  motion  for  an  atmospheric  reentry  vehicle  about  a  spher¬ 
ical,  rotating  Earth,  as  described  by  Vinh,13  with  coefficient  of  lift  (Cl)  and  coefficient  of 
drag  ( Cd )  are  used  to  describe  the  reentry  trajectory.  As  a  note,  the  subscript  d  denotes 
dimensional,  all  angles  remain  in  radians,  and  thrust  is  assumed  zero  as  vehicle  considered 
in  this  work  is  a  glide  vehicle.  Additionally,  terms  denoted  by  a  ’hat’,  x,  denote  a  nondi- 
mensional  representation  of  a  state  variable.  For  reference,  the  angles  used  in  the  equations 
of  motion  are  displayed  in  Figure  3. 
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Figure  3.  Spherical  Coordinates  for  Atmospheric  Flight:  longitude  6 ,  latitude  </>,  heading  angle  ip,  flight  path 
angle  7,  and  bank  angle  a 

First,  the  dimensional  kinematic  and  dynamic  equations  of  motion  are  presented  in  equa¬ 
tion  2  where  u>e  is  the  angular  velocity  of  earth’s  rotation. 

fd  =  Vd  sin  7 

■  Vd  cos  7  cos  0 

I'd  COS  0 

■  Vd  cos  7  sin  0 


Vd  = 


- - —  ()d  sin  7  +  TdOj'i  cos  0  [cos  0  sin  7  —  sin  0  sin  0  cos  7] 

md 

1  \Ld  ____  _  ____  ,  1 


7  =  —  —  cos  a  —  gd  cos  7  H - cos  7 

Vd  [  md  rd 

+2Vrfa;ecos0cos0  +  r^o;2  cos0  [cos0cos7  —  sin  0  sin  0  sin  7] 

0  =  —  — — - —  cos  7  cos0  tan0  +  2Vdue  [sin0  cos  0tan7  —  sin  < 

Vd[md  cos  J  rd 


sin  0  cos  0  cos  0 


For  the  purpose  of  optimization,  the  equations  in  2  are  normalised  so  that  a  change  in  one 
unit  of  each  state  element  is  approximately  of  equal  significance.14  These  equations  can  be 
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non-dimensionalized  and  simplified  using  the  variables  in  Eq.  (3). 13,15 
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(3) 


The  control  variables  are  bank  angle  (a)  and  the  angle  of  attack  (a),  which  ultimately 
regulates  the  lift  and  drag  terms  in  equation  2.  The  resulting  non-dimensional  equation  of 
motion  are: 
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(4) 


The  atmospheric  model  used  in  this  work  was  based  on  the  U.S.  Standard  Atmosphere 
in  1976,16,1'  with  an  extended  altitude  range  of  1000km.  Additionally,  a  spherical  earth  is 
considered  in  this  work  and  thus  results  in  a  Newtonian,  or  inverse  square  law,  gravity  field. 
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IV.  Results 


The  equations  presented  in  the  previous  section  are  used  to  replicate  generic  hyper¬ 
sonic  vehicle  simulation  results  derived  from  Program  to  Optimize  Simulated  Trajectories 
(POST).18,19  The  constant  speed/altitude  2-D  optimal  solution  was  previously  presented10 
and  derived  using  the  discontinuous  costate  dynamic  optimization  method.20  The  same 
technique  will  now  be  applied  to  the  3-dimensional  spherical  Earth  problem. 

This  section  compares  the  results  for  several  simulated  3-dimensional  optimal  trajectories 
as  calculated  by  both  POST  and  GPOCS.  A  set  of  desired  mission  profiles  are  established 
to  use  as  a  common  comparison  between  the  two  methods.  The  missions  range  in  constraint 
complexity,  beginning  with  a  benign  constraint  case  and  ending  with  a  heavily  constrained 
mission  profile.  The  desired  trajectories  consist  several  constraints,  which  include  control, 
heating,  altitude,  waypoint,  and  flight  restriction  zone  constraints.  The  results  of  this  anal¬ 
ysis  compare  and  contrast  the  effectiveness  of  the  optimization  methods  in  both  POST  and 
GPOCS  to  arrive  at  an  optimal  solution  for  each  of  the  constrained  mission  profiles.  The 
results  in  this  section  provide  insight  into  the  advantages  and  disadvantages  in  using  dif¬ 
ferent  optimization  methods  for  a  range  of  mission  complexities.  Three  comparisons  make 
the  overall  evaluation  of  the  two  methods;  the  first  to  compare  the  underlying  system  dy¬ 
namics  of  the  two  methods,  the  second  to  evaluate  a  simple  max  range  comparison,  and  the 
final  to  compare  a  heavily  constrained  trajectory.  The  evaluation  will  begin  with  the  first 
comparison  to  evaluate  the  system  dynamics. 

IV. A.  Dynamics  Comparison 

The  first  comparison  focuses  on  evaluating  the  basic  dynamics,  atmospheric  model,  and 
gravity  model  between  POST  and  GPOCS.  This  comparison  is  setup  to  ensure  subsequent 
comparisons  are  made  using  identical  system  dynamics,  atmospheric  model,  and  gravity 
model.  This  comparison  also  includes  an  evaluation  between  the  implicit  and  explicit  inte¬ 
gration  methods  used  in  GPOCS  and  POST,  respectively.  This  is  to  ensure  the  integration 
methods  deliver  identical  results  for  an  identical  problem  formulation.  As  stated  earlier,  the 
extended  US76  atmosphere  and  Newtonian  gravity  field  is  used  in  both  GPOCS  and  POST. 
To  ensure  identical  dynamics  setups  a  simple  open  loop,  non  optimised  trajectory  was  inte¬ 
grated  forward  in  time  using  a  single  fixed  angle  of  attack  and  bank  angle  in  both  POST  and 
GPOCS.  This  test  resulted  in  identical  trajectories  between  POST  and  GPOCS.  As  a  note, 
this  comparison  used  the  Runge-Kutta  fourth  order  integration  method19  for  both  POST 
and  the  dynamics  used  in  GPOCS.  GPOCS,  however,  used  an  implicit  integration  method, 
as  discussed  earlier.  Therefore,  before  any  comparison  in  optimization  was  to  be  done,  a 
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comparison  of  the  separate  integration  methods  was  be  performed.  The  initial  conditions 
and  results  of  such  comparison  are  presented  in  Figures  4  and  5.  These  results  indicated 
the  implicit  and  explicit  integration  methods  of  POST  and  GPOCS,  respectively,  produced 
identical  trajectory  results.  As  a  note,  this  test  used  open  loop,  non-optimised  trajectories, 
with  fixed  angles  of  attack  of  10  degrees  and  a  bank  angle  of  60  degree. 
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Figure  4.  Initial  Conditions  of  Integration  Method  Comparison 
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Latitude  versus  Longitude 


Figure  5.  Trajectory  Comparison  of  Integration  Method  Comparison 


IV.B.  Basic  Optimization  Comparison 

The  second  set  of  comparisons  involved  an  evaluation  of  the  optimization  methods  in 
both  GPOCS  and  POST.  For  this,  a  simple  maximum  downrange  trajectory  optimization 
was  performed.  The  first  trajectory  to  compare  with  the  pseudo-spectral  method  was  a 
straight  gliding  trajectory  meeting  the  requirements  of  the  DARPA  FALCON  program21  (ref: 
9000nm  glide  range)  for  a  20001b  high  lift  hypersonic  glide  vehicle.  The  description  of  this 
vehicle  is  outlined  in  Appendix  A.  For  the  POST  trajectory,  three  angles  of  attack  were  used 
to  control  the  glide  profile  and  allow  for  constraints  on  altitude  if  desired:  the  first  AOA 
during  reentry  determines  the  depth  into  the  atmosphere  the  hypersonic  vehicle  reenters 
before  ‘skipping’  up  into  a  glide,  the  second  AOA  determines  how  high  the  vehicle  ‘skips’, 
and  the  third  is  held  constant  for  the  remainder  of  the  trajectory  causing  a  ‘phugoid’  or 
‘skipping’  trajectory.  The  initial  conditions  (or  reentry  conditions)  were  found  given  the 
vehicle  weight  and  range  requirement.  The  resulting  reentry  velocity  and  flight  path  angle 
were  consistent  with  the  original  DARPA  FALCON  concept:  launching  a  20001b  HTV-1  from 
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a  SpaceX  Falcon  1  launch  vehicle.  The  three  initial  guesses  for  angles  of  attack  were  based 
on  experience.  Once  the  initial  conditions  and  guesses  were  set,  non-optimized  trajectories 
using  constant  alpha  and  then  constant  bank  angle  were  run  to  validate  the  equations  of 
motion  in  the  pseudo-spectral  model.  To  ensure  simplicity  in  this  evaluation  the  only  criteria 
that  was  evaluated  was  imply  the  cost  of  the  final  longitude  value,  with  the  only  constraints 
being  a  minimum  altitude  value  of  27.4  kilometers,  or  90k  feet,  and  a  final  time  of  3063 
seconds,  that  being  the  run  time  of  the  POST  trajectory.  The  maximum  time  was  chosen 
to  more  closely  align  the  GPOCS  results  to  that  of  the  POST  run.  The  initial  conditions, 
altitude,  attained  down  range  values,  and  angle  of  attack  comparison  plots  are  presented  in 
Figures  6,7,8,  and  9. 
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Figure  6.  Initial  Conditions  of  the  First  Optimization  Comparison 
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Altitude 


Figure  7.  Altitude  Comparison  of  the  First  Optimization  Comparison 


Inertial  Coordinates 
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Figure  8.  Range  Comparison  of  the  First  Optimization  Comparison 
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Angle  of  Attack 


Figure  9.  Angle  of  Attack  Comparison  of  the  First  Optimization  Comparison 


The  greatest  distinction  revealed  from  this  comparison  was  the  flexability  in  available 
control  afforded  to  the  user  in  GPOCS,  as  seen  in  Figure  9  .  Herein,  the  user  defines  a 
control  input  vector  containing  the  value  of  each  control  at  every  node  along  the  trajectory. 
This  is  contrasted  with  POST,  where  the  user  can  only  specify  discrete  control  inputs  which 
are  then  optimised  to  attain  the  maximum  down  range  value.  Consequently,  finer  control 
variability  in  GPOCS  has  resulted  in  a  greater  downrange  value,  as  seen  in  Figure  8.  The 
next  comparison  will  present  a  more  constrained  problem  to  compare  and  contrast. 

IV. C.  Complex  Optimization  Comparison 

In  this  comparison,  POST  was  set  up  to  fly  a  hypersonic  glide  trajectory  following  a 
ballistic  reentry  with  fixed  reentry  conditions.  Here,  a  problem  was  set  up  to  launch  a 
hypersonic  glide  vehicle  from  Cape  Canaveral  to  the  Persian  Gulf,  flying  over  the  Mediter¬ 
ranean  Sea  with  minimal  overflight  of  land.  This  was  done  by  specifying  a  waypoint  at 
the  Strait  of  Gibraltar  and  two  no-fly  zones,  one  in  North  Africa  and  one  in  South-central 
Europe.  POST  simulated  the  glide  from  reentry  so  an  arbitrary  reentry  point  (see  output 
hies)  was  determined  as  an  initial  condition.  The  reentry  velocity  and  flight  path  angle 
were  kept  the  same  as  the  prior  case,  but  reentry  azimuth  was  allowed  to  vary  assuming 
a  booster  would  tailor  its  profile  to  meet  the  desired  end  conditions.  The  trajectory  was 
optimized  for  maximum  final  altitude,  terminating  at  a  set  velocity.  This  forced  POST  to 
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converge  to  a  solution  that  would  conserve  the  most  energy,  i.e.  not  waste  energy  through 
excessive  maneuvering.  However,  significant  maneuvering  was  required  to  hit  the  waypoint 
and  avoid  the  no-fly  zones.  Considering  that  reentry  azimuth  could  vary,  this  required  three 
bank  maneuvers  for  each  geographic  constraint  and  a  fourth  bank  maneuver  to  straighten 
the  glide  before  the  terminal  condition.  Bank  angles,  bank  duration,  and  time  between  turns 
were  set  up  as  independent  variables  to  go  along  with  the  three  angles  of  attack  from  before. 
Dependent  variables  were  reentry  angle,  and  the  geographic  constraints  of  the  waypoint  and 
no-fly  zones.  These  constraints  were  defined  by  using  ‘tracking  stations’  in  POST  and  apply¬ 
ing  limits  to  ‘slant  range’  either  forcing  the  trajectory  towards  or  away  from  the  particular 
station. 

To  set  up  a  problem  such  as  this  in  POST,  one  cannot  apply  all  the  constraints  at  once 
and  expect  POST  to  come  to  a  solution.  The  problem  must  be  built  up  as  constraints  are 
applied  one  at  a  time.  Using  the  trajectory  from  the  pseudo-spectral  validation  and  adding 
four  banking  maneuvers,  one  could  quickly  converge  to  a  solution  that  when  projected  on 
a  Mercator  Projected  map,  looks  close  to  the  desired  solution  .  However,  the  guesses  for 
angles  of  attack,  bank  angles,  and  durations  were  all  based  on  experience  and  did  not  taking 
into  account  the  geographic  constraints.  This  ‘first-cut’  is  useful  however  for  generating 
new  independent  variables  to  be  used  once  these  constraints  are  applied.  Additionally  by 
outputting  the  slant  ranges  for  each  of  the  three  tracking  stations,  even  though  they  are  not 
yet  used  for  constraints,  one  can  make  changes  in  the  independent  variables  to  observe  how 
the  slant  ranges  change  to  better  determine  the  first  guess  once  the  constraints  are  applied. 

Even  though  the  above  described  trajectory  is  close  to  the  solution,  POST  cannot  yet 
handle  all  three  geographic  constraints  being  applied  at  once.  Therefore  the  problem  added 
the  complexity  of  the  waypoint  located  at  the  Strait  of  Gibraltar  and  the  user  could  continue 
to  monitor  the  slant  ranges  to  the  other  two  tracking  stations.  This  process  continued,  adding 
one  geographic  constraint  at  a  time  until  all  constraints  were  applied  and  the  problem  could 
be  optimized.  Even  then,  POST  was  unable  to  converge  on  a  ‘best’  solution  although  all  the 
constraints  were  met. 

Once  the  first  trajectory  to  the  Persian  Gulf  was  found  using  POST  each  geographic 
constraint  was  applied  one  at  a  time.  The  geographic  constraints  were  based  on  tracking 
stations  with  slant  ranges  used  to  either  keep  the  trajectory  near  or  far  from  the  stations. 
Because  of  the  narrow  width  of  the  Strait  of  Gibraltar,  the  tolerance  for  the  waypoint  was 
first  a  quarter  nautical  mile,  but  this  prevented  POST  from  converging  on  a  solution.  The 
tolerance  was  then  relaxed  to  one  nautical  mile  and  POST  was  then  able  to  converge  on 
a  solution.  The  no-fly  zones  were  defined  to  avoid  overflight  of  land  in  North  Africa  and 
South-central  Europe,  so  a  minimum  distance  to  each  tracking  station  was  defined.  Here, 
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the  tolerance  on  each  slant  range  had  to  be  relaxed  to  five  significant  digits  ,  as  the  slant 
range  is  measured  in  feet.  This  translated  to  approximately  one  a  two-thirds  nautical  miles. 

The  first  trajectory  to  the  Persian  Gulf  without  any  geographic  constraints  took  a  CPU 
time  of  .4  seconds  .  Applying  the  first  waypoint  increased  the  CPU  time  to  .5  sec  .  The 
first  no-fly  zone  made  the  CPU  time  at  .6  sec.  This  case  reached  the  iteration  limit  without 
finding  a  solution,  ffowever,  by  noting  the  errors  on  the  constraints,  the  user  was  able 
to  modify  the  independent  variables  to  make  a  reasonable  first  guess  with  all  geographic 
constraints  applied.  Even  so,  this  case  also  ran  to  the  iteration  limit  without  converging  on 
a  solution.  In  this  case  CPU  time  was  still  .6  seconds.  In  all,  this  process  took  an  experienced 
POST  user  (with  experience  modeling  hypersonic  glide  trajectories,  enabling  educated  first 
guesses)  a  day  and  a  half,  running  POST  approximately  30  times  with  up  to  50  iterations 
each. 

In  comparison,  the  method  to  generate  the  mission  described  above  using  GPOCS  was 
significantly  easier.  The  overall  trajectory  as  compared  to  POST  is  presented  in  Figure  10. 


Figure  10.  Initial  Conditions  of  the  Final  Optimization  Comparison 


Upon  the  user  defining  the  waypoint,  being  the  end  of  the  first  phase,  and  the  no  fly 
zones  as  the  path  constraints  the  trajectory  is  ready  to  be  run.  Simply  put,  no  iterative 
method  is  necessary  to  generate  the  trajectory  as  in  the  case  of  POST.  The  results  of  the 
trajectory  comparison  is  seen  in  Figures  11  through  15 
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Velocity,  km/sec  Latitude,  deg  Altitude, 


Altitude 


Figure  11.  Altitude  Comparison  of  the  Final  Optimization  Comparison 


Latitude  versus  Longitude 


Figure  12.  Latitude  and  Longitude  of  the  Final  Optimization  Comparison 


Velocity 


Figure  13.  Velocity  of  the  Final  Optimization  Comparison 
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Total  Heating 


Figure  14.  Total  Heat  of  the  Final  Optimization  Comparison 


Stagnation  Heat  Flux 


Figure  15.  Peak  Heating  of  the  Final  Optimization  Comparison 


V.  Conclusion 

This  paper  presents  an  example  of  employing  the  Gauss  Pseudospectral  Optimization 
method  in  optimal  trajectory  generation  of  a  hypersonic  glide  vehicle.  A  brief  description 
of  the  implicit  optimization  method  used  in  the  Pseudospectral  method  is  discussed  and 
contrasted  with  the  method  used  in  POST.  A  series  of  comparison  optimal  trajectories  are 
presented  and  reviewed.  The  comparisons  reveal  GPOCS  generated  trajectories  produce 
similar  results  to  those  of  POST.  However,  the  iterative  approach  commonly  used  in  POST 
is  not  necessary  to  successfully  generate  a  complex  trajectory  in  GPOCS. 
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A.  Hypersonic  Vehicle  Aerodynamic  Data 

The  following  is  taken  from18  and  is  used  to  create  a  Mach  independent  model  for  use 
this  research. 


Table  A.l.  High  Lift  Hypersonic  Glide  Vehicle  Aero  Data  Base 


Lift  to  Drag  Ratio  ( L/ D ) 

AOA 

Mach  3.5 

Mach  5 

Mach  8 

Mach  10 

Mach  15 

Mach  20 

Mach  23 

10° 

2.2000 

2.5000 

3.1000 

3.5000 

3.3846 

3.2692 

3.2000 

15° 

2.5000 

2.6616 

2.9846 

3.2000 

3.0846 

2.9692 

2.9000 

20° 

2.2000 

2.3616 

2.6846 

2.9000 

2.7846 

2.6692 

2.6000 

Coefficient  of  Lift  {Cl) 

AOA 

Mach  3.5 

Mach  5 

Mach  8 

Mach  10 

Mach  15 

Mach  20 

Mach  23 

10° 

0.4500 

0.4250 

0.4000 

0.3800 

0.3700 

0.3600 

0.3500 

15° 

0.7400 

0.7000 

0.6700 

0.6300 

0.6000 

0.5700 

0.5570 

20° 

1.0500 

1.0000 

0.9500 

0.9000 

0.8500 

0.8000 

0.7800 

Coefficient  of  Drag  {Co) 

AOA 

Mach  3.5 

Mach  5 

Mach  8 

Mach  10 

Mach  15 

Mach  20 

Mach  23 

10° 

0.2045 

0.1700 

0.1290 

0.1090 

0.1090 

0.1090 

0.1090 

15° 

0.2960 

0.2630 

0.2240 

0.1970 

0.1950 

0.1920 

0.1920 

20° 

0.4770 

0.4230 

0.3540 

0.3100 

0.3050 

0.3000 

0.3000 

HGV-H  Aero  Reference  Area  Sref  =  750  in2 
HGV-H  Mass  m  =  2000  lbs  =  907.186  kg  (assuming  g  =  32.174  ft/s2) 
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B.  Gauss  Pseudospectral  Method  Exposition 


For  completeness,  an  introduction  to  the  Gauss  Pseudospectral  Method  (GPM),  which 
is  the  basis  for  the  software  package  used  in  this  research,  is  presented  below.  The  GPM  is 
an  orthogonal  collocation  method  where  the  collocation  points  are  the  Legendre- Gauss  (LG) 
points.  The  description  herein  is  a  compilation  from1, 22-26  which  is  based  on.27,28  Additional 
implementation  methods  are  in.26,29-31 


B.I.  Continuous  Bolza  Problem 


The  dynamic  optimization  problem  is  restated  here  with  the  transformation  of  the  inde¬ 
pendent  variable  t: 

t  r  -  t  r»  f  r 

(B.I) 


tf  -t0  tf  +  t0 

t  =  - - - T  H - - — 


The  optimal  control  problem  is  to  determine  the  state,  x(r)  G  Mn,  control,  u(r)  G  Mm,  initial 
time,  to,  and  final  time,  tf,  that  minimizes  the  cost  functional: 


J  =  0(x(-l),  t0,  x(l ),tf)  +  tf  2  to  J  L(x(r),  u(r),  r;  t0,  tf)dr  (B.2) 

subject  to  the  constraints 

“j~T  =  ~  2  — /(x(r),u(T),r;G,t/)  (B.3) 

^(x(-l),t0,x(l),t/)  =  0  (B.4) 

C(x(r),  u(t),  r;  tQ,  tf)  <  0  (B.5) 


Herein,  the  optimal  control  problem  of  Eqs.  (B.2-B.5)  is  called  the  continuous  Bolza  problem. 


B.II.  Gauss  Pseudospectral  Discretization  of  Continuous  Bolza  Problem 

The  direct  approach  to  solving  the  continuous  Bolza  optimal  control  problem  of  Sec.  B.I 
is  to  discretize  and  transcribe  Eqs.  (B.2-B.5)  to  a  nonlinear  programming  problem  (NLP). 
The  Gauss  pseudospectral  method,  like  Legendre  and  Chebyshev  methods,  is  based  on  ap¬ 
proximating  the  state  and  control  trajectories  using  interpolating  polynomials.  The  state  is 
approximated  using  a  basis  of  N  +  1  Lagrange  interpolating  polynomials,  C. 

N 

*M«X(T)  =  £x(7i)  £(t)  (B.6) 

i- 0 
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where  £;(r)(i  =  0, . . . ,  N)  are  defined  as 


N 

aw  =  n 

j=0,j^i 


T  -  T-j 


n  -  r,- 


(B.7) 


Additionally,  the  control  is  approximated  using  a  basis  of  iV  Lagrange  interpolating  polyno¬ 
mials  £*(r),  (i  =  1, . . . ,  N)  as 


iV 


u(r)  «  U(r)  =  ^U(r,;)£* 


i= 1 


where 


N 


r*M  =  11 


r  -  Ti 


,  Tj  —  T: 
j=l,j^i  3 


It  can  be  seen  from  Eqs.  (B.7)  and  (B.9)  that  Ci{r)(i  —  0, . . . ,  N)  and  £*(r)(f 
satisfy  the  properties 


(B.8) 

(B.9) 
!,•••,  N) 


A(Tj)  = 


1  ,  i  =  j 
0  ,  i^j 


and  £*(tj)  = 


1  ,  i  =  j 

0  ,  i^j 


(B.10) 


Differentiating  the  expression  in  Eq.  (B.6)  produces 


N 


x(t)  «  X(t)  =  ^2  x(a)A(t) 


(B.ll) 


i=0 


The  derivative  of  each  Lagrange  polynomial  at  the  LG  points  can  be  represented  in  a  differ¬ 
ential  approximation  matrix,  D  G  W.NxN+l.  The  elements  of  the  differential  approximation 
matrix  are  determined  offline  as  follows: 


N 


N 


n  (r^_ri) 


N 

'  II 

j=0,j^i 


(B.12) 


Ti  Tj ) 
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where  k  —  1, . . . ,  N  and  i  =  0, . . . ,  N.  The  dynamic  constraint  is  transcribed  into  algebraic 
constraints  via  the  differential  approximation  matrix  as  follows: 

N  t  -t 

Y  Dki*i  -  -^y -/(Xfc,  Ufc,  Tk-  tQ,tf)  =  0  (k  =  1, . . . ,  N)  (B.13) 

4= 0 

where  Xk  =  X(rk)  G  and  U&  =  U(rjfc)  G  (k  =  1, . . . ,  N).  Note  that  the  dynamic 
constraint  is  collocated  only  at  the  LG  points  and  not  at  the  boundary  points  (this  form  of 
collocation  differs  from  other  well  known  pseudospectral  methods31,32).  Additional  variables 
in  the  discretization  are  defined  as  follows:  X0  =  X(— 1),  and  X/-  =  X(l),  where  Xf  is 
defined  in  terms  of  Xfc,  (k  —  0, . . . ,  N)  and  U (rk)(k  —  1, . . . ,  N)  via  the  Gauss  quadrature.33 

t  -t  N 

X/  =  X o  +  wkf(Kk,  Ufc,  Tfc,;  t0,tf)  (B.14) 

k= 1 

The  continuous  cost  function  of  Eq.  (B.2)  is  approximated  using  a  Gauss  quadrature33  as 

,  ,  * 

J  =  0(XO,  to,  Xf,  tf )  +  — -  "Y;  wkL(Xk,  Ufc,  Tfc.;  to,  tf )  :  (B.15) 

k= 1 

where  wk  are  the  Gauss  weights.  The  boundary  constraint  of  Eq.  (B.4)  is  expressed  as 

^(Xo,to,Xf,tf)  —  0  (B.16) 

Furthermore,  the  path  constraint  of  Eq.  (B.5)  is  evaluated  at  the  LG  points  as 


C(Xk,XJk,Tk;t0,tf)  <  0  (k  =  l,...,N) 


(B.17) 


The  cost  function  of  Eq.  (B.15)  and  the  algebraic  constraints  of  Eqs.  (B.13),  (B.14),  (B.16), 
and  (B.17)  define  an  NLP  whose  solution  is  an  approximate  solution  to  the  continuous  Bolza 
problem.  Finally,  it  is  noted  that  the  above  discretization  can  be  employed  in  multiple- 
phase  problems  by  transcribing  the  problem  in  each  phase  using  the  above  discretization 
and  connecting  the  phases  by  linkage  constraints: 


p(«)(x(pf)(i/),^);q(pf),x(pi)(i0),iJ,S);q(pS))  =  0,  (phpu  G  [1, . . . ,  P] ,  s  =  1, 


5  ^PJ 


(B.18) 

where  x^(f)  G  Wnp,u^p\t)  G  Mmp,q^  G  and  t  G  M  are,  respectively,  the  state,  control, 
static  parameters,  and  time  in  phase  p  G  [1, . . . ,  P],  Lp  is  the  number  of  phases  to  be  linked, 
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pf  G  [1, . . . ,  P] ,  (s  —  1, ,  Lp)  are  the  “left”  phase  numbers,  and  p®  G  [1, . . . ,  P] ,  (s  — 
1, . ... ,  Lp)  are  the  “right”  phase  numbers,  where  the  sum  of  the  cost  for  each  phase: 

p 

J  =  Y l  J(p)  (B.19) 

p= i 


becomes  the  final  cost  to  be  minimized. 

B.III.  KKT  Conditions  of  the  NLP 

The  hrst-order  optimality  conditions  (i.  e.  ,  the  Karush- Kuhn- Tucker  (KKT)  condi¬ 
tions)  of  the  NLP  can  be  obtained  using  the  augmented  cost  function  or  Lagrangian.  The 
augmented  cost  function  is  formed  using  the  Lagrange  multipliers  Afc  G  Mn,  jlk  G  Mc, 
k  =  1, .  . . ,  N,  Af  G  Mn,  and  u  G  W  as 

t  t  N 

Ja  =  0(Xo,  to,  X/,  tf)  +  — -  22  WkL(X.k,  Ufc,  Tk ;  to,  t/)  —  UT,lp(K 0,  to,  X/,  tf ) 

fc=l 

Af  N  /  N 

—  22  Ufc,  rk ;  tG,  tf)  —  22  A[  I  ^  DkiXi  —  ,f  -  —  /(Xfc,  Ufc,  rfc;  t0,  tf) 

k= 1  fc=l  \i=l 

-  Al  ( Xf  -  X0  -  ^  wfc/(Xfc,  Ufc,  rfc;  to,  t/) 

V  fc=i 

(B.20) 

The  KKT  conditions  are  found  by  setting  equal  to  zero  the  derivatives  of  the  Lagrangian 
with  respect  to  X0,  Xfc,  Xy,  Ufc,  Afc,  jlk,  AF,  u,  t0,  tf.  The  solution  to  the  NLP  of  Sec.  B.ll 
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must  satisfy  the  following  KKT  conditions: 


N 


E  X‘D 


tf  ~t 0, 
ki  0  Jk 


i= 0 


N  (!T 

E  t  +  + 

i=l  V 


tf-to  (  dLk 


5Xz 


A-^  + A[ 

wk 


t 

fc7V+l 


5/fe  2  jll  dCk 


5Xfc  tf  — 10  wk  5Xfc 


o  =  ^±  + 

5Uz.  + 


A[  ,  At \  dfk  2  ftdCk 
wk  FJ  d\Jk  tf-t0wkdUk 

0(XO,  t0,  Xy,  tf)  =  0 
t  t1  dd)  ~~rT'  dib 

\T0=-—  +  VT 


<9X 


o 


ax, 


o 


~  t-i  dd>  _r  ah' 
Ah  =  — z/  — 


X 


JV 


/  */ 
N 


tf  —  to  ^Hk  1  TT  ~T 

-^^E“,‘757  +  9E^  =  75:-i/  757 


fc=l 


0 


k= 1 


0 


0 


tf  —  to  <97/*,  ,  1  u 

2^  Wk7F7  +  9  2^  Wfci^ 


fc=i 


dtf 

7  fc=i 

Ufc  <0 

Jljk  =  0,  when  Cjk  <  0 
fljk  <  0,  when  Cjfc  =  0 


50  _r  50 

a7  +  " 


Ap  —  A0  + 


tf  -  t0 


AT 

E 

k= 1 


Wfc 


x/  =  x0 


5Lfc 

5Xi 


tf  ~ 


AT 


^  ^  IT  kfk 


k= 1 

'  Ar 

—  +  aJ 

wk 


5/fe 


2  fil  8C ) 


+ 


tf  to  Wh  d~X.k 


(B.21) 

(B.22) 

(B.23) 

(B.24) 

(B.25) 

(B.26) 

(B.27) 

(B.28) 

(B.29) 

(B.30) 

(B.31) 

(B.32) 

(B.33) 


where  the  shorthand  notation  Lk  =  L(Xk,  Ufc,  rfc;  t0,  t/),  fk  =  /(Xfc,  Ufc,  rfc;  t0,  t/),  = 

htfc(Xfc,  Afe,  }lk,  Ufc,  Tfc;  t0,  t/),  and  Cjk  =  Cj(Xk,  Ufe,  rfc;  t0,  t/)  is  used.  The  augmented  Hamil¬ 
tonian,  Hk,  is  defined  as 


Hk  =  Tfc  + 


fk  ~ 


2 

- Ofc 

tf  to  tCfc 


(B.34) 
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and  An  is  defined  as 


dtp  _T  dip 


9X0  '  '  <9X0  (B’35) 

Note  that  the  sign  of  the  jij.  term  is  opposite  to  the  pr  term;  thus,  the  software  output 
requires  a  sign  change  to  match  the  analytical  results  presented  prior  to  this  Appendix. 


B.IV.  First-Order  Optimality  Conditions  of  Continuous  Bolza  Problem 

The  indirect  approach  to  solving  the  continuous  Bolza  problem  of  Eqs.  (B.2-B.5)  in 
Sec.  B.l  is  to  apply  the  calculus  of  variations  and  Pontryagin’s  Maximum  Principle34  to 
obtain  the  first-order  necessary  conditions  for  optimality.35  These  variational  conditions  are 
typically  derived  using  the  augmented  Hamiltonian,  H ,  defined  as 


H(x,  A,  n,  u,  r;  t0,tf)  =  L(x,  u,  r;  t0,tf)  +  Ar(r)/(x,  u,  r;  t0,tf)  -  fiT(r)C(x,  u,  r;  t0,tf) 

(B.36) 

where  A (r)  G  R"  is  the  costate  and  /t(r)  G  Rc  is  the  Lagrange  multiplier  associated  with  the 
path  constraint.  The  continuous-time  first-order  optimality  conditions  can  be  shown  to  be 


tf  ~  to  ,,  ,  ,  \  tf  —t o  OH 

s-—  /(X-U'T;  k'tf)  =  —-§i 


dA 

dr 


dx 

dr 
tf  -  t0 


dL 

dx 


X 


dx 


+  [I 


dC  A  tf  —  to  dH 

fa)  =  I 


dx 


dL  .  rr,  df  rj-\  dC 


_  dH 

du  du 

^(x(To  ),t0,x(Tf),tf)  =  0 


0  -  k  +  ATai 


A(r0)  = 


dp 


+  v 


dip 


dx(r0)  dx(r0 )  ’ 


Mh)  = 


dp 


v 


dp 


«w- 


dx(Tf)  dx(Tf) 

dp  T  dp 


I-Ij(t)  =  0,  when  Cj(x,  u,  r;  t0,  tf)  <  0,  j  =  1, . . . ,  c 
Ab(r)  <  0,  when  C,-(x,  u,  r;  t0,  tf)  =  0,  j  =  1, . . . ,  c 


(B.37a) 

(B.37b) 

(B.37c) 

(B.37d) 

(B.37e) 

(B.37f) 

(B.37g) 

(B.37h) 


where  v  G  R9  is  the  Lagrange  multiplier  associated  with  the  boundary  condition  p.  It  can 
be  shown  that  the  augmented  Hamiltonian  at  the  initial  and  final  times  can  be  written, 
respectively,  as 

HW  =  --ptl  £  |br  +  \  Mr  (B.38) 

^«/)  =  ^  £  f^dr  +  5  /_;  Mr  (B.39) 
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B.V.  Gauss  Pseudospectral  Discretization  Necessary  Conditions 

In  order  to  discretize  the  variational  conditions  of  Sec.  B.IV  using  the  Gauss  pseudospec¬ 
tral  discretization,  it  is  necessary  to  form  a  suitable  approximation  for  the  costates.  In  this 
method,  the  costate,  A (r),  is  approximated  as  follows: 

JV+l 

A(r)  «  A(r)  =  Y  A(T*)At(T)  (B-40) 

i= 1 


where  C\(r){i  —  1, . . . ,  N  +  1)  are  defined  as 


n+i 

dw=  n  (b-41) 

j=l,j^i  *  3 

It  is  emphasized  that  the  costate  approximation  is  different  from  the  state  approximation.  In 
particular,  the  basis  of  the  IV +  1  Lagrange  interpolating  polynomials  C\{r){i  =  1, . . . ,  IV  +  1) 
includes  the  costate  at  the  final  time  (as  opposed  to  the  initial  time  which  is  used  in  the 
state  approximation).  This  (non- intuitive)  costate  approximation  is  necessary  in  order  to 
provide  a  complete  mapping  between  the  KKT  conditions  and  the  variational  conditions. 

Using  the  costate  approximation  of  Eq.  (B.40),  the  first-order  necessary  conditions  of  the 
continuous  Bolza  problem  in  Eq.  (B.37)  are  discretized  as  follows.  First,  the  state  and  control 
are  approximated  using  Eq.  (B.6)  and  Eq.  (B.8),  respectively.  Next,  the  costate  is  approxi¬ 
mated  using  the  basis  of  IV  +  1  Lagrange  interpolating  polynomials  as  defined  in  Eq.  (B.40). 
The  continuous-time  first-order  optimality  conditions  of  Eq.  (B.37)  are  discretized  using  the 
variables  X0  =  X(— 1),  X^.  =  X(rk)  G  Rn,  and  Xj  =  X(l)  for  the  state;  U*,  =  U(rfc)  G  Rm 
for  the  control;  A0  =  A(— 1),  Ak  =  A [rk)  G  R",  and  AF  =  A(l)  for  the  costates;  and 
fik  =  n{rk)  G  Rc,  for  the  Lagrange  multiplier  associated  with  the  path  constraints  at  the  LG 
points  k  =  1 , ,N.  The  other  unknown  variables  in  the  problem  are  the  initial  and  final 
times,  to  £  R,  £/  G  R,  and  the  Lagrange  multiplier,  v  G  RT  The  total  number  of  variables 
is  then  given  as  (2 n  +  m  +  c)N  +  4n  +  q  +  2.  These  variables  are  used  to  discretize  the 
continuous  necessary  conditions  of  Eq.  (B.37)  via  the  Gauss  pseudospectral  discretization. 
Note  that  the  derivative  of  the  state  is  approximated  using  Lagrange  polynomials  based  on 
IV  +  1  points  consisting  of  the  N  LG  points  and  the  initial  time,  To,  while  the  derivative  of 
the  costate  is  approximated  using  Lagrange  polynomials  based  on  IV  +  1  points  consisting  of 
the  N  LG  points  and  the  final  time,  Tf  .  The  resulting  algebraic  equations  that  approximate 
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the  continuous  necessary  conditions  at  the  LG  points  are  given  as 
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(B.42) 

(B.43) 

(B.44) 

(B.45) 

(B.46) 

(B.47) 

(B.48) 

(B.49) 

(B.50) 

(B.51) 


for  k  =  1, . . . ,  N  and  j  =  1 , ,c.  The  final  two  equations  that  are  required  (in  order  to  link 
the  initial  and  final  state  and  costate,  respectively)  are 


Af  —  A0  + 


f  _  +  N 

Xf  =  X0  +  ^Y^Y^wkfk 

k= 1 

tf  —  to  (  dLk  T  dfk  T  dCk  \ 

—  X^{-mr  K*mCt+ 


(B.52) 

(B.53) 


The  total  number  of  equations  in  the  set  of  discrete  necessary  conditions  of  Eqs.  (B.42-B.53) 
is  (2 n  +  m  +  c)N  +  4n  +  q  +  2  (the  same  number  of  unknown  variables).  Solving  these 
nonlinear  algebraic  equations  would  be  an  indirect  solution  to  the  optimal  control  problem. 


B.VI.  Costate  Estimate 

One  of  the  key  features  of  the  Gauss  pseudospectral  method  is  the  ability  to  map  the  KKT 
multipliers  of  the  NLP  to  the  costates  of  the  continuous-time  optimal  control  problem.  In 
particular,  using  the  results  of  Sections  B.III  and  B.V,  a  costate  estimate  for  the  continuous 
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Bolza  problem  can  be  obtained  at  the  Legendre- Gauss  points  and  the  boundary  points. 
This  costate  estimate  is  taken  from27  and  is  summarized  below  via  the  Gauss  Pseudospectral 
Costate  Mapping  Theorem-. 

Theorem  B.l.  Gauss  Pseudospectral  Costate  Mapping  Theorem:  The  Karush- 
Kuhn-Tucker  (KKT)  conditions  of  the  NLP  are  exactly  equivalent  to  the  discretized  form  of 
the  continuous  first-order  necessary  conditions  of  the  continuous  Bolza  problem  when  using 
the  Gauss  pseudospectral  discretization.  Furthermore,  a  costate  estimate  at  the  initial  time, 
final  time,  and  the  Legendre- Gauss  points  can  be  found  from  the  KKT  multipliers,  Ak,  jlk, 
AF,  and  v , 

Ak  =  — 1  +  AF,  gk  =  - — — ,  v  =  v,  A  (t0)  =  A0,  A  (tf)  =  AF  (B.54) 

wk  tf-t0wk 

LIsing  the  substitutions  of  Eq.  (B.54),  it  can  be  seen  that  Eqs.  (B.21-B.33)  are  exactly  the 
same  as  Eqs.  (B.42-B.53). 

B.VII.  Computation  of  Boundary  Controls 

It  is  seen  in  the  GPM  that  the  control  is  discretized  only  at  the  LG  points  and  is  not 
discretized  at  either  the  initial  or  the  terminal  point.  Consequently,  the  solution  of  the  NLP 
defined  by  Eqs.  (B.13-B.17)  does  not  produce  values  of  the  controls  at  the  boundaries.  The 
ability  to  obtain  accurate  initial  and  terminal  controls  can  be  important  in  many  applications, 
particularly  in  guidance  where  real-time  computation  of  the  initial  control  is  of  vital  interest. 

At  first  glance,  it  may  seem  that  the  lack  of  control  information  at  the  boundaries  can 
be  overcome  simply  via  extrapolation  of  the  control  at  the  LG  points.  However,  multiple 
reasons  exist  as  to  why  this  is  not  the  best  approach.  First,  no  particular  functional  form  for 
the  control  is  assumed  in  the  GPM  discretization.  As  a  result,  the  best  function  to  use  for 
extrapolation  is  ambiguous.  Second,  any  reasonable  extrapolation  of  the  control  (e.g.,  linear, 
quadratic,  cubic,  or  spline)  may  violate  a  path  constraint  which,  in  general,  will  render  the 
extrapolated  control  infeasible.  Third,  even  if  the  extrapolated  control  is  feasible,  it  will 
not  satisfy  the  required  optimality  conditions  at  the  boundaries  (i.  e.  the  control  will  be 
suboptimal  with  respect  to  the  NLP).  Consequently,  it  is  both  practical  and  most  rigorous 
to  develop  a  systematic  procedure  to  compute  the  boundary  controls.  The  primal  (state) 
and  dual  (costate)  solutions  of  the  NLP  arising  from  the  Gauss  pseudospectral  method  are 
used  to  compute  the  boundary  controls.25 

Computation  of  the  initial  control  is  done  first  since  the  approach  for  computing  the  ter¬ 
minal  control  is  identical.  First,  recalling  the  augmented  Hamiltonian,  H ,  for  the  continnous- 
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time  optimal  control  problem  in  Eq.  (B.36)  is 


tf(x,u,A,/x)  =  L  +  \T f  -  jiTC  (B.55) 

where  shorthand  notation  is  used.  Recall  from  the  principle  of  Pontryagin,  at  every  instant 
of  time  the  optimal  control  is  the  control  u*(r)  G  U  that  satisfies  the  condition 

LT(x*,  u* ,  A*,//)  <  LT(x*,  u,  A*, /V)  (B.56) 

where  U  is  the  feasible  control  set.  Consequently,  for  a  given  instant  of  time  r  where  x*(t), 
A*(t),  and  /i*(r)  are  known,  Eq.  (B.56)  is  a  constrained  optimization  problem  in  u(r)  G  Mm. 
In  order  to  solve  this  constrained  optimization  problem  at  the  initial  time,  it  is  necessary  to 
know  x*(to),  A*(to),  and  /x*(ro). 

Consider  the  information  that  can  be  obtained  by  solving  the  NLP  associated  with  the 
GPM.  In  particular,  the  primal  solution  to  the  NLP  produces  X(to)  while  the  dual  solution 
to  the  NLP  can  be  manipulated  algebraically  to  obtain  the  initial  costate,  A(r0).  However, 
because  the  NLP  does  not  evaluate  the  path  constraint  at  the  boundaries,  there  is  no  asso¬ 
ciated  Lagrange  multiplier  This  apparent  impediment  can  be  overcome  by  applying 

the  minimum  principle  in  a  manner  somewhat  different  from  that  given  in  Eq.  (B.56).  In 
particular,  suppose  we  let  H  be  the  Hamiltonian  (not  the  augmented  Hamiltonian),  where 
7i  is  defined  as 

W(x,  u,  A)  =  L  +  \Tf  (B.57) 

It  is  seen  in  Eq.  (B.57)  that  the  term  involving  the  path  constraint  is  not  included.  The 
path  constraint  is  instead  incorporated  into  the  feasible  control  set.  In  particular,  suppose 
we  let  Vo  be 

Vo  =W  p|  Co  (B.58) 

where  Vo  is  the  intersection  of  the  original  set  of  feasible  controls  at  time  To,  denoted  7/,  with 
the  set  of  all  controls  at  time  r0  that  satisfy  the  inequality  constraint  of  Eq.  (B.17),  denoted 
Co-  Then,  using  the  values  of  X(r0)  and  A(t0),  the  following  modified  optimization  problem 
in  m  variables  U(to)  G  Mm  can  be  solved  to  obtain  the  initial  control,  U(to): 

min  7d(X(r0),  U(r0),  A(r0),  r0;  t0,  tf)  (B.59) 

U(r0)GV0 

It  is  noted  that,  because  Vo  is  restricted  by  the  inequality  path  constraint  at  r0,  the  solution 
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of  U(t0)  is  equivalent  to  the  solution  of  the  following  problem: 


min  W(X(r0),U(ro),A(ro),r0;to,t/) 

U(r0)GW  (B.60) 

subject  to  C(X(t0),  U(t0),  t0;  t0,tf)  <  0 

Interestingly,  if  the  constraint  is  active,  then  the  initial  path  constraint  multiplier,  /i*(to),  will 
also  be  determined  by  the  minimization  problem  of  Eq.  (B.60).  Finally,  similar  to  the  initial 
time,  the  control  at  the  terminal  time,  U (ry),  can  be  obtained  by  solving  the  minimization 
problem  of  Eq.  (B.60)  at  r  =  77,  i.  e. 

m\n7  U(rj),  A(t/),  rf;  t0,  tf) 

uh/)ew  (B.61) 

subject  to  C(X.(Tf),XJ(Tf),Tf;to,tf)  <  0 

The  Gauss  Pseudospectral  Optimal  Control  Software  (GPOCS)  is  a  software  program  written 
in  MATLAB3  for  solving  multiple-phase  optimal  control  problems,  which  implements  the 
algorithm  as  described  above. 


aMATLAB  is  a  registered  trademark  of  The  Mathworks,  Inc.,  3  Apple  Hill  Drive,  Natick,  MA 
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